home *** CD-ROM | disk | FTP | other *** search
/ Aminet 35 / Aminet 35 (2000)(Schatztruhe)[!][Feb 2000].iso / Aminet / game / shoot / ADescentSrc.lha / descent / fix / fix.c next >
C/C++ Source or Header  |  1998-03-29  |  12KB  |  541 lines
  1. /*
  2. THE COMPUTER CODE CONTAINED HEREIN IS THE SOLE PROPERTY OF PARALLAX
  3. SOFTWARE CORPORATION ("PARALLAX").  PARALLAX, IN DISTRIBUTING THE CODE TO
  4. END-USERS, AND SUBJECT TO ALL OF THE TERMS AND CONDITIONS HEREIN, GRANTS A
  5. ROYALTY-FREE, PERPETUAL LICENSE TO SUCH END-USERS FOR USE BY SUCH END-USERS
  6. IN USING, DISPLAYING,  AND CREATING DERIVATIVE WORKS THEREOF, SO LONG AS
  7. SUCH USE, DISPLAY OR CREATION IS FOR NON-COMMERCIAL, ROYALTY OR REVENUE
  8. FREE PURPOSES.  IN NO EVENT SHALL THE END-USER USE THE COMPUTER CODE
  9. CONTAINED HEREIN FOR REVENUE-BEARING PURPOSES.  THE END-USER UNDERSTANDS
  10. AND AGREES TO THE TERMS HEREIN AND ACCEPTS THE SAME BY USE OF THIS FILE.  
  11. COPYRIGHT 1993-1998 PARALLAX SOFTWARE CORPORATION.  ALL RIGHTS RESERVED.
  12. */
  13. /*
  14.  * $Source: /usr/CVS/descent/fix/fix.c,v $
  15.  * $Revision: 1.7 $
  16.  * $Author: tfrieden $
  17.  * $Date: 1998/03/28 23:06:13 $
  18.  * 
  19.  * C version of fixed point library
  20.  * 
  21.  * $Log: fix.c,v $
  22.  * Revision 1.7  1998/03/28 23:06:13  tfrieden
  23.  * fixmulaccum moved to assembler, now function pointer
  24.  *
  25.  * Revision 1.6  1998/03/27 00:54:37  tfrieden
  26.  * Got rid of this monster overkill fixdivquadlong
  27.  *
  28.  * Revision 1.5  1998/03/25 21:59:02  tfrieden
  29.  * Removed the long int stuff for now
  30.  *
  31.  * Revision 1.4  1998/03/22 19:21:19  tfrieden
  32.  * added new command line arguments for fix point math
  33.  *
  34.  * Revision 1.3  1998/03/22 15:21:53  tfrieden
  35.  * function pointers finally work
  36.  *
  37.  * Revision 1.2  1998/03/22 01:49:47  tfrieden
  38.  * Inserted function pointer for fixmul
  39.  *
  40.  * Revision 1.1.1.1  1998/03/03 15:11:50  nobody
  41.  * reimport after crash from backup
  42.  *
  43.  * Revision 1.2  1998/02/28 01:14:24  tfrieden
  44.  * Don`t really know what changed
  45.  *
  46.  * Revision 1.1.1.1  1998/02/13  20:20:12  hfrieden
  47.  * Initial Import
  48.  */
  49.  
  50.  
  51. #pragma off (unreferenced)
  52. static char rcsid[] = "$Id: fix.c,v 1.7 1998/03/28 23:06:13 tfrieden Exp $";
  53. #pragma on (unreferenced)
  54.  
  55. #include <stdlib.h>
  56. #include <math.h>
  57.  
  58. #include "error.h"
  59. #include "fix.h"
  60.  
  61. extern ubyte guess_table[];
  62. extern short sincos_table[];
  63. extern ushort asin_table[];
  64. extern ushort acos_table[];
  65. extern fix isqrt_guess_table[];
  66.  
  67. // stuff for fixpoint arithemtic pointers
  68. // standard stuff
  69. fix fixmul_double(fix a __asm("d0"), fix b __asm("d1"));
  70. fix fixdiv_double(fix a __asm("d0"), fix b __asm("d1"));
  71. fix fixmuldiv_double(fix a __asm("d0"), fix b __asm("d1"), fix c __asm("d2"));
  72. long fixdivquadlong_double(ulong nl __asm("d0"),ulong nh __asm("d1"),ulong d __asm("d2"));
  73. // integer versions using long instructions
  74. extern fix fixmul_int(fix a __asm("d0"), fix b __asm("d1"));
  75. extern fix fixdiv_int(fix a __asm("d0"), fix b __asm("d1"));
  76. extern long fixdivquadlong_int(ulong nl __asm("d0"),ulong nh __asm("d1"),ulong d __asm("d2"));
  77. void fixmulaccum_int(quad *q __asm("a0"),fix a __asm("d0"),fix b __asm("d1"));
  78. // obsolete
  79. // extern fix fixmul_short(fix a __asm("d0"), fix b __asm("d1"));
  80. // FPU versions
  81. extern fix fixmul_fpu(fix a __asm("d0"), fix b __asm("d1"));
  82. extern fix fixdiv_fpu(fix a __asm("d0"), fix b __asm("d1"));
  83. extern fix fixmuldiv_fpu(fix a __asm("d0"), fix b __asm("d1"), fix c __asm("d2"));
  84. extern long fixdivquadlong_fpu(ulong nl __asm("d0"),ulong nh __asm("d1"),ulong d __asm("d2"));
  85. void fixmulaccum_fpu(quad *q __asm("a0"),fix a __asm("d0"),fix b __asm("d1"));
  86. // Function pointers for the above functions
  87. fix (*fixmul_func)(fix a __asm("d0"), fix b __asm("d1"));
  88. fix (*fixdiv_func)(fix a __asm("d0"), fix b __asm("d1"));
  89. fix (*fixmuldiv_func)(fix a __asm("d0"), fix b __asm("d1"), fix c __asm("d2"));
  90. long (*fixdivquadlong_func)(ulong nl __asm("d0"),ulong nh __asm("d1"),ulong d __asm("d2"));
  91. void (*fixmulaccum_func)(quad *q __asm("a0"),fix a __asm("d0"),fix b __asm("d1"));
  92.  
  93.  
  94. #undef DEBUG_PROFILE
  95. #ifdef DEBUG_PROFILE
  96. extern fix profile_fa2_time, profile_fd_time, profile_fm_time, profile_fmd_time;
  97. extern int profile_fa2_called, profile_fd_called, profile_fm_called, profile_fmd_called;
  98. #endif
  99.  
  100. void fix_init(void)
  101. {
  102.  
  103.     fixmul_func         = fixmul_int;
  104.     fixdiv_func         = fixdiv_int;
  105.     fixmuldiv_func      = fixmuldiv_double;
  106.     fixdivquadlong_func = fixdivquadlong_int;
  107.     fixmulaccum_func    = fixmulaccum_int;
  108.  
  109.     if (FindArg("-fpu")) {
  110.         fixmul_func     = fixmul_fpu;
  111.         fixdiv_func     = fixdiv_fpu;
  112.         fixmuldiv_func  = fixmuldiv_fpu;
  113.         fixdivquadlong_func = fixdivquadlong_fpu;
  114.         fixmulaccum_func    = fixmulaccum_fpu;
  115.     }
  116.  
  117. }
  118.  
  119. //negate a quad
  120. void fixquadnegate(quad *q)
  121. {
  122.     q->low  = 0 - q->low;
  123.     q->high = 0 - q->high - (q->low != 0);
  124. }
  125.  
  126. //multiply two ints & add 64-bit result to 64-bit sum
  127. //void fixmulaccum_old(quad *q __asm("a0"),fix a __asm("d0"),fix b __asm("d1"))
  128. void fixmulaccum_old(quad *q, fix a,fix b)
  129. {
  130.     ulong aa,bb;
  131.     ulong ah,al,bh,bl;
  132.     ulong t,c=0,old;
  133.     int neg;
  134.  
  135.     neg = ((a^b) < 0);
  136.  
  137.     aa = labs(a); bb = labs(b);
  138.  
  139.     ah = aa>>16;  al = aa&0xffff;
  140.     bh = bb>>16;  bl = bb&0xffff;
  141.  
  142.     t = ah*bl + bh*al;
  143.  
  144.     if (neg)
  145.         fixquadnegate(q);
  146.  
  147.     old = q->low;
  148.     q->low += al*bl;
  149.     if (q->low < old) q->high++;
  150.     
  151.     old = q->low;
  152.     q->low += (t<<16);
  153.     if (q->low < old) q->high++;
  154.     
  155.     q->high += ah*bh + (t>>16) + c;
  156.     
  157.     if (neg)
  158.         fixquadnegate(q);
  159.  
  160. }
  161.  
  162. //extract a fix from a quad product
  163. fix fixquadadjust(quad *q)
  164. {
  165.     return (q->high<<16) + (q->low>>16);
  166. }
  167.  
  168.  
  169. #define EPSILON (F1_0/100)
  170.  
  171. fix fixmul_double(fix a __asm("d0"), fix b __asm("d1"))
  172. {
  173.     return (fix)(((double)a*(double)b)/65536.0);
  174. }
  175.  
  176. fix fixdiv_double(fix a __asm("d0"), fix b __asm("d1"))
  177. {
  178.     return (fix)(((double)a * 65536.0) / (double)b);
  179. }
  180.  
  181. fix fixmuldiv_double(fix a __asm("d0"), fix b __asm("d1"), fix c __asm("d2"))
  182. {
  183.     double d;
  184.  
  185.     d = (double)a * (double) b;
  186.     return (fix)(d / (double) c);
  187.  
  188. }
  189.  
  190. //given cos & sin of an angle, return that angle.
  191. //parms need not be normalized, that is, the ratio of the parms cos/sin must
  192. //equal the ratio of the actual cos & sin for the result angle, but the parms 
  193. //need not be the actual cos & sin.  
  194. //NOTE: this is different from the standard C atan2, since it is left-handed.
  195.  
  196. fixang fix_atan2(fix cos,fix sin)
  197. {
  198.     double d, dsin, dcos;
  199.     fix m,t;
  200.  
  201.     //Assert(!(cos==0 && sin==0));
  202.  
  203.     //find smaller of two
  204.  
  205.     dsin = (double)sin;
  206.     dcos = (double)cos;
  207.     d = sqrt((dsin * dsin) + (dcos * dcos));
  208.  
  209.     if (d==0.0) {
  210.         return 0;
  211.     }
  212.  
  213.     if (labs(sin) < labs(cos)) {                //sin is smaller, use arcsin
  214.         t = fix_asin((fix)((dsin / d) * 65536.0));
  215.         if (cos<0)
  216.             t = 0x8000 - t;
  217.         return t;
  218.     }
  219.     else {
  220.         t = fix_acos((fix)((dcos / d) * 65536.0));
  221.         if (sin<0)
  222.             t = -t;
  223.         return t;
  224.     }
  225.  
  226. }
  227.  
  228. //divide a quad by a fix, returning a fix
  229. long fixdivquadlong_double(ulong nl __asm("d0"),ulong nh __asm("d1"),ulong d __asm("d2"))
  230. {
  231.     int i;
  232.     ulong tmp0;
  233.     ubyte tmp1;
  234.     ulong r;
  235.     ubyte T,Q,M;
  236.  
  237.     r = 0;
  238.  
  239.     Q = ((nh&0x80000000)!=0);
  240.     M = ((d&0x80000000)!=0);
  241.     T = (M!=Q);
  242.  
  243.     if (M == 0)
  244.     {
  245.         for (i=0; i<32; i++ )   {
  246.     
  247.             r <<= 1;
  248.             r |= T;
  249.             T = ((nl&0x80000000L)!=0);
  250.             nl <<= 1;
  251.     
  252.             switch( Q ) {
  253.         
  254.             case 0:
  255.                 Q = (unsigned char)((0x80000000L & nh) != 0 );
  256.                 nh = (nh << 1) | (unsigned long)T;
  257.  
  258.                 tmp0 = nh;
  259.                 nh -= d;
  260.                 tmp1 = (nh>tmp0);
  261.                 if (Q == 0)
  262.                     Q = tmp1;
  263.                 else
  264.                     Q = (unsigned char)(tmp1 == 0);
  265.                 break;
  266.             case 1:
  267.                 Q = (unsigned char)((0x80000000L & nh) != 0 );
  268.                 nh = (nh << 1) | (unsigned long)T;
  269.  
  270.                 tmp0 = nh;
  271.                 nh += d;
  272.                 tmp1 = (nh<tmp0);
  273.                 if (Q == 0)
  274.                     Q = tmp1;
  275.                 else
  276.                     Q = (unsigned char)(tmp1 == 0);
  277.                 break;
  278.             }
  279.             T = (Q==M);
  280.         }
  281.     }
  282.     else
  283.     {
  284.         for (i=0; i<32; i++ )   {
  285.     
  286.             r <<= 1;
  287.             r |= T;
  288.             T = ((nl&0x80000000L)!=0);
  289.             nl <<= 1;
  290.     
  291.             switch( Q ) {
  292.         
  293.             case 0:
  294.                 Q = (unsigned char)((0x80000000L & nh) != 0 );
  295.                 nh = (nh << 1) | (unsigned long)T;
  296.  
  297.                 tmp0 = nh;
  298.                 nh += d;
  299.                 tmp1 = (nh<tmp0);
  300.                 if (Q == 1)
  301.                     Q = tmp1;
  302.                 else
  303.                     Q = (unsigned char)(tmp1 == 0);
  304.                 break;
  305.             case 1: 
  306.                 Q = (unsigned char)((0x80000000L & nh) != 0 );
  307.                 nh = (nh << 1) | (unsigned long)T;
  308.  
  309.                 tmp0 = nh;
  310.                 nh = nh - d;
  311.                 tmp1 = (nh>tmp0);
  312.                 if (Q == 1)
  313.                     Q = tmp1;
  314.                 else
  315.                     Q = (unsigned char)(tmp1 == 0);
  316.                 break;
  317.             }
  318.             T = (Q==M);
  319.         }
  320.     }
  321.  
  322.     r = (r << 1) | T;
  323.  
  324.     return r;
  325. }
  326.  
  327. ulong quad_sqrt(long low,long high)
  328. {
  329.     long cnt,r,old_r,t;
  330.     quad tq;
  331.  
  332.     if (high<0)
  333.         return 0;
  334.  
  335.     if (high==0 && low>=0)
  336.         return long_sqrt(low);
  337.  
  338.     if (high & 0xff000000)
  339.         cnt=12+16;
  340.     else if (high & 0xff0000)
  341.         cnt=8+16;
  342.     else if (high & 0xff00)
  343.         cnt=4+16;
  344.     else
  345.         cnt=0+16;
  346.     
  347.     r = guess_table[(high>>cnt)&0xff]<<cnt;
  348.  
  349.     //quad loop usually executed 4 times
  350.  
  351.     r = (fixdivquadlong(low,high,r)+r)/2;
  352.     r = (fixdivquadlong(low,high,r)+r)/2;
  353.     r = (fixdivquadlong(low,high,r)+r)/2;
  354.  
  355.     do {
  356.  
  357.         old_r = r;
  358.         t = fixdivquadlong(low,high,r);
  359.  
  360.         if (t==r)   //got it!
  361.             return r;
  362.  
  363.         r = (t+r)/2;
  364.  
  365.     } while (!(r==t || r==old_r));
  366.  
  367.     t = fixdivquadlong(low,high,r);
  368.     tq.low=tq.high;
  369.     fixmulaccum(&tq,r,t);
  370.     if (tq.low!=low || tq.high!=high)
  371.         r++;
  372.  
  373.     return r;
  374. }
  375.  
  376. //computes the square root of a long, returning a short
  377. ushort long_sqrt(long a)
  378. {
  379.     int cnt,r,old_r,t;
  380.  
  381.     if (a<=0)
  382.         return 0;
  383.  
  384.     if (a & 0xff000000)
  385.         cnt=12;
  386.     else if (a & 0xff0000)
  387.         cnt=8;
  388.     else if (a & 0xff00)
  389.         cnt=4;
  390.     else
  391.         cnt=0;
  392.     
  393.     r = guess_table[(a>>cnt)&0xff]<<cnt;
  394.  
  395.     //the loop nearly always executes 3 times, so we'll unroll it 2 times and
  396.     //not do any checking until after the third time.  By my calcutations, the
  397.     //loop is executed 2 times in 99.97% of cases, 3 times in 93.65% of cases, 
  398.     //four times in 16.18% of cases, and five times in 0.44% of cases.  It never
  399.     //executes more than five times.  By timing, I determined that is is faster
  400.     //to always execute three times and not check for termination the first two
  401.     //times through.  This means that in 93.65% of cases, we save 6 cmp/jcc pairs,
  402.     //and in 6.35% of cases we do an extra divide.  In real life, these numbers
  403.     //might not be the same.
  404.  
  405.     r = ((a/r)+r)/2;
  406.     r = ((a/r)+r)/2;
  407.  
  408.     do {
  409.  
  410.         old_r = r;
  411.         t = a/r;
  412.  
  413.         if (t==r)   //got it!
  414.             return r;
  415.  
  416.         r = (t+r)/2;
  417.  
  418.     } while (!(r==t || r==old_r));
  419.  
  420.     if (a % r)
  421.         r++;
  422.  
  423.     return r;
  424. }
  425.  
  426. //computes the square root of a fix, returning a fix
  427. fix fix_sqrt(fix a)
  428. {
  429.     return ((fix) long_sqrt(a)) << 8;
  430. }
  431.  
  432.  
  433. //compute sine and cosine of an angle, filling in the variables
  434. //either of the pointers can be NULL
  435. //with interpolation
  436. void fix_sincos(fix a,fix *s,fix *c)
  437. {
  438.     int i,f;
  439.     fix ss,cc;
  440.  
  441.     i = (a>>8)&0xff;
  442.     f = a&0xff;
  443.  
  444.     ss = sincos_table[i];
  445.     if (s) *s = (ss + (((sincos_table[i+1] - ss) * f)>>8))<<2;
  446.  
  447.     cc = sincos_table[i+64];
  448.     if (c) *c = (cc + (((sincos_table[i+64+1] - cc) * f)>>8))<<2;
  449. }
  450.  
  451. //compute sine and cosine of an angle, filling in the variables
  452. //either of the pointers can be NULL
  453. //no interpolation
  454. void fix_fastsincos(fix a,fix *s,fix *c)
  455. {
  456.     int i;
  457.  
  458.     i = (a>>8)&0xff;
  459.  
  460.     if (s) *s = sincos_table[i] << 2;
  461.     if (c) *c = sincos_table[i+64] << 2;
  462. }
  463.  
  464. //compute inverse sine
  465. fixang fix_asin(fix v)
  466. {
  467.     fix vv;
  468.     int i,f,aa;
  469.  
  470.     vv = labs(v);
  471.  
  472.     if (vv >= f1_0)     //check for out of range
  473.         return 0x4000;
  474.  
  475.     i = (vv>>8)&0xff;
  476.     f = vv&0xff;
  477.  
  478.     aa = asin_table[i];
  479.     aa = aa + (((asin_table[i+1] - aa) * f)>>8);
  480.  
  481.     if (v < 0)
  482.         aa = -aa;
  483.  
  484.     return aa;
  485. }
  486.  
  487. //compute inverse cosine
  488. fixang fix_acos(fix v)
  489. {
  490.     fix vv;
  491.     int i,f,aa;
  492.  
  493.     vv = labs(v);
  494.  
  495.     if (vv >= f1_0)     //check for out of range
  496.         return 0;
  497.  
  498.     i = (vv>>8)&0xff;
  499.     f = vv&0xff;
  500.  
  501.     aa = acos_table[i];
  502.     aa = aa + (((acos_table[i+1] - aa) * f)>>8);
  503.  
  504.     if (v < 0)
  505.         aa = 0x8000 - aa;
  506.  
  507.     return aa;
  508. }
  509.  
  510. #define TABLE_SIZE 1024
  511.  
  512. //for passed value a, returns 1/sqrt(a) 
  513. fix fix_isqrt( fix a )
  514. {
  515.     int i, b = a;
  516.     int cnt = 0;
  517.     int r;
  518.  
  519.     if ( a == 0 ) return 0;
  520.  
  521.     while( b >= TABLE_SIZE )    {
  522.         b >>= 1;
  523.         cnt++;
  524.     }
  525.  
  526.     //printf( "Count = %d (%d>>%d)\n", cnt, b, (cnt+1)/2 );
  527.     r = isqrt_guess_table[b] >> ((cnt+1)/2);
  528.  
  529.     //printf( "Initial r = %d\n", r );
  530.  
  531.     for (i=0; i<3; i++ )    {
  532.         int old_r = r;
  533.         r = fixmul( ( (3*65536) - fixmul(fixmul(r,r),a) ), r) / 2;
  534.         //printf( "r %d  = %d\n", i, r );
  535.         if ( old_r >= r ) return (r+old_r)/2;
  536.     }
  537.  
  538.     return r;   
  539. }
  540.  
  541.